建立微分方程模型一般可分以下三步：
- （1）根据实际要求确定研究的量（自变量、未知函数、必要的参数等），并确定坐标系。
- （2）找出这些量所满足的基本规律。
- （3）运用这些规律列出方程和定解条件。 下面通过实例介绍常用的利用微分方程建立数学模型的方法之微元分析法
微元分析法¶
该方法的基本思想是通过分析研究对象的有关变量在一个很短时间内的变化规律，寻找一些微元之间的关系式。
问题¶
有高为1m的半球形容器，水从它的底部小孔流出。小孔横截面积为1cm$^2$。开始时容器内盛满了水，求水从小孔流出过程中容器里水面的高度 （水面与孔口中心的距离）随时间 变化的规律。
建模¶
如图所示, 以底部中心为坐标原点, 垂直向上为坐标轴的正向建立坐标系。 由水力学知, 水从孔口流出的流量 $Q$ 为 “通过孔口横截面的水的体积 $V$ 对时间 $t$ 的变化率", 满足 $$ Q=\frac{d V}{d t}=0.62 S \sqrt{2 g h}, $$ 式中, $0.62$ 为流量系数; $g$ 为重力加速度 (取 $9.8 \mathrm{~m} / \mathrm{s}^{2}$ ), $S$ 为孔口横截面积 (单位: $\left.\mathrm{m}^{2}\right), \boldsymbol{h}$ 为 $t$ 时刻水面高度 (单位: $\mathrm{cm}$ )。
求解¶
当 $S=1 \mathrm{~cm}^{2}=0.0001 \mathrm{~m}^{2}$ 时, $$ d V=0.000062 \sqrt{2 g h} d t . $$ 在微小时间间隔 $[t, t+d t]$ 内, 水面高度由 $h$ 降到 $h+d h$ (这里 $d h<0$ ), 容 器中水的体积改变量近似为 $$ d V=-\pi r^{2} d h, $$ 式中, $r$ 为 $t$ 时刻的水面半径, 右端置负号是由于 $d h<0$, 而 $d V>0$; 这里 $$ r^{2}=\left[1^{2}-(1-h)^{2}\right]=2 h-h^{2} . $$
由上式, 得 $0.000062 \sqrt{2 g h} d t=\pi\left(h^{2}-2 h\right) d h$. 再考虑到初始条件, 得到如下的微分方程模型 $$ \left\{\begin{array}{l} \frac{d t}{d h}=\frac{10000 \pi}{0.62 \sqrt{2 g}}\left(h^{\frac{3}{2}}-2 h^{\frac{1}{2}}\right), \\ t(1)=0, \end{array}\right. $$ 利用分离变量法, 可以求得微分方程的解为 $$ t(h)=-15260.5042 h^{\frac{3}{2}}+4578.1513 h^{\frac{5}{2}}+10682.3530 . $$
代码
import sympy as sp
sp.var('h') #定义符号变量
sp.var('t', cls=sp.Function) #定义符号函数
g = 9.8
eq = t(h).diff(h) -10000*sp.pi/0.62/sp.sqrt(2*g)*(h**(3/2)-2*h**(1/2)) #定义方程
t = sp.dsolve(eq, ics={t(1): 0}) #求微分方程的符号解
t = sp.simplify(t)
print(t.args[1].n(9))
参考文献
- 司守奎，孙玺菁 Python数学实验与建模